#include "tdma.h"
#include <assert.h>

////////////////////////////////////////////////
///////////////PRIVATE MEMBER VARIABLES/////////
////////////////////////////////////////////////
void tdma::Algorithm( )
{
   
   double tempc, tempd;
   double xTemp;
   //UPPER-DIAGONAL MATRIX
   
   for ( int i=1;i<n;i++ )
   {
      //c values      
      tempc = gsl_vector_get( c, i ) - 
            gsl_vector_get( b, i )*gsl_vector_get( c, i-1 )/
            gsl_vector_get( d, i-1 );
      gsl_vector_set( c, i, tempc );      
      // d values
      tempd = gsl_vector_get( d, i ) - 
            gsl_vector_get( b, i )*gsl_vector_get( a, i-1 )/
            gsl_vector_get( d, i-1 );
      gsl_vector_set( d, i, tempd );
   }
   //BACK SUBSTITUTION
   double numerator;
   double denominator;
   //n-1 th variable
   xTemp = gsl_vector_get( c, n-1 )/gsl_vector_get( d, n-1 );
   //cout<<"Xtemp :"<<xTemp<<endl;
   gsl_vector_set( x, n-1, xTemp );
   //all other variables
   for ( int i=n-2;i>=0;i-- )
   {
      numerator = ( gsl_vector_get( c, i ) - gsl_vector_get( a, i )*gsl_vector_get( x, i+1 ) );
      denominator = gsl_vector_get( d, i );
      xTemp = numerator/denominator;
      //cout<<"Xtemp :"<<xTemp<<endl;
      gsl_vector_set( x, i, xTemp ); 
   }
}

////////////////////////////////////////////////
///////////////PUBLIC MEMBER VARIABLES//////////
////////////////////////////////////////////////
tdma::tdma( )
{
}


void tdma::Solve( gsl_vector* bInp, gsl_vector* dInp, gsl_vector* aInp, 
      gsl_vector* cInp )
{
   n = bInp->size;
   if ( DBG_TDMA ) cout<<"Size is :"<<bInp->size<<endl;
   //set bool variables
   sqmatset = false;
   //make sure vector sizes are the same on all vectors   
   assert( bInp->size = dInp->size );
   assert( bInp->size = aInp->size );
   assert( bInp->size = cInp->size );
   
   //assign variable ti inputs
   b = bInp;
   d = dInp;
   a = aInp;
   c = cInp;
   x = gsl_vector_calloc( n );
   //solve the system
   Algorithm( );
}

void tdma::Delete( )
{
   //cout<<"cleaning up tdma :"<<endl;
   gsl_vector_free( x );

}

tdma::~tdma( )
{
   cout<<"cleaning up tdma :"<<endl;
   gsl_vector_free( b );
   gsl_vector_free( d );
   gsl_vector_free( a );
   gsl_vector_free( c );
   if ( sqmatset == true )
      gsl_matrix_free( SqMat );
   gsl_vector_free( x );
}

/*tdma::tdma( gsl_matrix* sqmat, gsl_vector* cInp )
{
   SqMat = gsl_matrix_calloc( sqmat->size1, sqmat->size2 );
   sqmatset = true;
   //make sure it is a square matrix
   assert( SqMat->size1 = SqMat->size2 );
   n = SqMat->size1;
   
   //assign variable ti inputs
   //c = gsl_vector_calloc( n );
   c = cInp;
   SqMat = sqmat;
   //allocate memory
   //b = gsl_vector_calloc( n );
   //d = gsl_vector_calloc( n );
   //a = gsl_vector_calloc( n );
   
   //get elements from the square matrix and put them into corresponding vectors
   double bTemp, dTemp, aTemp;   
   for ( int i=0;i<n;i++ )
   {
      //traverse the diagonal
      dTemp = gsl_matrix_get( SqMat, i, i );
      gsl_vector_set( d, i, dTemp );
      if ( i != 0 )
         
      {
         bTemp = gsl_matrix_get( SqMat, i, i-1 );
         gsl_vector_set( b, i, bTemp );
      }
      if ( i!= n-1 )
      {
         aTemp = gsl_matrix_get( SqMat, i, i+1 );
         gsl_vector_set( a, i, aTemp );
         
      }
   }
   //solve the system
   Algorithm( );
}*/

gsl_vector* tdma::GetX( )
{
   
   return x;
}


